/*************
author:: libo
license:: MIT
*****************/

#ifndef LB_Math_H_
#define LB_Math_H_
#define factorial LB_factorial
#define SIGN(x) ((x>0)?1:((x<0)?-1:0))
#ifndef SAFE_SQRT
#define SAFE_SQRT(x) (x>=0?sqrt(x):0)
#endif
#ifndef SAFE_FREE
#define SAFE_FREE(x) if(x!=NULL){free(x);x=NULL;}
#endif

// #ifndef DOUBLEZERO
// #define DOUBLEZERO(x)  (fabs(x)<=x)
// #endif

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<string.h>

#ifdef __cplusplus
extern "C" {
#endif

static inline  int double_equal_zero(double x)
{
    return (fabs(x)<=1e-12);
}

int factorial(int);
int * pailie_map_i_to_reverse_order(int);
int* map_reverse_order_to_pailie(int*,int);

// 在n个元素中选取m个，总共有多少种方案
// re_len是方案的个数
//
int* compute_combination(int m,int n,int* re,int* re_len);

static inline double *out_product(double *a,double *b)
{
    double* re=(double*)malloc(sizeof(double)*3);
    re[0]=a[1]*b[2]-a[2]*b[1];
    re[1]=-a[0]*b[2]+a[2]*b[0];
    re[2]=a[0]*b[1]-a[1]*b[0];
    return re;
}

static inline void out_productn(double *a,double *b,double *re)
{
    re[0]=a[1]*b[2]-a[2]*b[1];
    re[1]=-a[0]*b[2]+a[2]*b[0];
    re[2]=a[0]*b[1]-a[1]*b[0];
    //return re;
}
double inner_product(double* a,double* b,int len );

int normalize(double *a,int len);

// double lb_norm(double* p,int n);

//  result 可以等于a1 可以等于a2
//
static inline void array_plus(double *a1,double* a2,double* result,int len)
{
    for(int i=0;i<len;i++)
    {
        result[i]=a1[i]+a2[i];
    }
}
//  result 可以等于a1 可以等于a2

static inline void array_minus(double *a1,double* a2,double* result,int len)
{
    for(int i=0;i<len;i++)
    {
        result[i]=a1[i]-a2[i];
    }
}
//
static inline void array_mult_scala(double* a1,double fac,double* result,int len)
{
    for(int i=0;i<len;i++)
    {
        result[i]=a1[i]*fac;
    }
}

static double lb_norm(double* p,int n)
{
    double re=0;
    for(int i=0;i<n;i++)
    {
        re+=p[i]*p[i];
    }
    return SAFE_SQRT(re);
}

static double distance_of_two_points(double * a,double* b,int n)
{
    double re=0;
    for(int i=0;i<n;i++)
    {
        re+=(a[i]-b[i])*(a[i]-b[i]);
    }
    return SAFE_SQRT(re);
}
// double distance_of_two_points(double * a,double* b,int n);
// tent映射，均匀分布
static inline double tent_map(double x)
{
    return x<0.5?(2.0*x):(2.0*(1-x));
}
// sine map
static inline double sin_map(double x)
{
    return sin(M_PI*x);
}



#ifdef __cplusplus
}
#endif
#endif
#undef factorial
